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Abstract 

We derive a collisionless kinetic theory for an ensemble of molecules undergoing nonholo- 
nomic rolling dynamics. We demonstrate that the existence of nonholonomic constraints 
leads to problems in generalizing the standard methods of statistical physics. In particular, 
we show that even though the energy of the system is conserved, and the system is closed in 
the thermodynamic sense, some fundamental features of statistical physics such as invariant 
measure do not hold for such nonholonomic systems. Nevertheless, we are able to construct 
a consistent kinetic theory using Hamilton's variational principle in Lagrangian variables, by 
regarding the kinetic solution as being concentrated on the constraint distribution. A cold 
fluid closure for the kinetic system is also presented, along with a particular class of exact 
solutions of the kinetic equations. 
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1 Introduction 

Constrained dynamical systems play an important role in modern mechanics and statistical physics. 
The constraints applied to the system can be separated into two broad categories - holonomic 
and nonholonomic. Holonomic constraints restrict the particle motion to lie a certain surface 
in the configuration space. Nonholonomic constraints are then defined as any constraint that 
cannot be reduced to motion on a particular surface in the configuration space. There are some 
classical examples of nonholonomic systems with the constraints that are linear in velocities. These 
systems usually (but not necessarily) come from perfect friction limitation, such as rolling particles 
[1] (Chaplygin's ball) or, more broadly, a particular connection between several components of 
velocities, such as Chaplygin's sleigh [2]. One may also see [3, 4, 5] for recent discussions of these 
type of problems. The Lagrange-D'Alembert principle [G] is usually used to treat the dynamics 
of such systems. We refer the reader to the book of Bloch [6] for a discussion of nonholonomic 
dynamics, as well as a review of recent literature and methods. If constraints are not linear in 
velocities, such as isokinetic models, typically. Gauss' minimal force principle is used to derive the 
equations of motion. 

In keeping with the spirit of regular statistical mechanics, one would like to develop a kinetic 
theory for large number of coupled nonholonomic particles, akin to the Vlasov or Boltzmann theory 
of interacting gas particles. However, in general, the presence of nonholonomic constraints destroys 
the Hamiltonian structure for the dynamics of individual particles. There are exceptions when the 
Hamiltonian structure of the dynamics can be restored, but these are special cases [7, 8]. Without 
an underlying Hamiltonian structure, the development of the kinetic theory for nonholonomically 
constrained systems is a formidable challenge. 

There certainly have been substantial developments in the study of the statistical physics of 
systems with nonholonomic constraints -the isokinetic models - which enforce the constant tem- 
perature condition for molecular particle simulations. The isokinetic restriction, which is quadratic 
in the particle velocities, cannot be solved by either the classical Lagrange-D'Alembert method or 
its generalizations such as the Hamilton- Pontryagin method [9]. Instead, the methods of minimal 
constraints due to Gauss has been used to describe the system; see [10, 11, 12] for for some recent 
progress and a review of the literature. A short discussion of this progress is warranted here. 

If, in the absence of constraints, the microscopic particle motion is described in phase space by 
the equation z = X{z), then the corresponding transport equation for the distribution function 
f{z, t) is taken to be of the form 



If it is assumed that the vector field X{z) is Hamiltonian with z = [p, q) and X = {Hp, —Hg), in 
which case it has zero divergence. The standard methods of statistical mechanics then apply, in 
particular, the conservation law of entropy S = j f log / holds. 

In case when vector field is not Hamiltonian, the situation is more complex. As was first realized 
in [13] (without proof), set in differential-geometric context in [14, 15] and further extended in 
[16, 17, 18, 19], the non-Hamiltonian vector fields lead, generally, to curved geometry of phase 
space. It was shown that a more advantageous, and geometrically correct, version of transport 
equation is obtained by introducing the metric y/g = exp J dwzX{z) dt into (1) as follows: 




(1) 



-(^f)+diVz{X{z)^f)=0 



(2) 
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One can then prove that the generahzed entropy S = J ^/g f log / is conserved. These gen- 
eral considerations were then successfully applied to the systems with constraints that break the 
Hamiltonian structure, in particular, to the non-holonomic isokinetic constrain enforcing constant 
temperature: 

-mi±^ = Kq = const . (3) 

i 

While this theory is general and is applicable to both linear and nonlinear constraints, a successful 
application of this theory hinges on the computation of the volume element ^Jg[z). Unfortunately, 
as we show below in Remark 4.6, this approach is not applicable for the case of interacting non- 
holonomic particles considered here. First, the explicit computation of the volume element is not 
possible, and second, and more important, there are persistent fundamental difficulties with proper 
definition of the divergences of vector field X for our case. Thus, unfortunately, we were not able 
to define a conservative entropy-like quantity, because the usual definition of entropy produces 
a divergent integral. This is perhaps because every particle in the ensemble is nonholonomically 
constrained to roll individually rather than having a single constraint for the whole system, such 
as (3). 

In this paper, we develop a nonholonomic kinetic theory for the particular case of interacting 
rolling particles whose centers of mass are offset from their geometrical centers and whose inertia 
tensors are not proportional to unity. Only rolling nonholonomic constraint is applied to the 
molecules, similar to recent work [ ] which treated this system as a model for investigating the 
properties of molecular monolayers. In that work, the particles had the same mechanical properties 
as the water molecules and thus their inertia tensors did not satisfy the symmetry requirements for 
the Chaplygin's ball. The theory we have developed here may also be viewed as an augmentation 
of the recently developed stochastic nonhonomic dynamics [ ]. The main contributions of this 
paper to the literature can be divided into three main topics. 

First main theme is the derivation of the equations of motion , and Sections 2,3 and 4 are 
devoted to that topic. We will emphasize that in principle, we could have derived the equations 
of motion and corresponding kinetic equation using the Gauss' method of minimal constraint. In 
that case, we would need to utilize recently developed geometric extension of this theory [22], valid 
for constraints in arbitrary spaces, and not only in M". However, we believe that such a derivation, 
although useful, would be exceedingly cumbersome; indeed, the Gauss' method of minimal con- 
straint force is usually used only when a small number of constraints is present in the system. In 
our case, it is important to emphasize that the rolling constraints are applied to every particle in 
the ensemble, so we have instead used the corresponding geometric generalization of previous ideas 
developed by [23] in the context of plasma physics with a Lagrange-d'Alembert's principle. We do 
not know of another work deriving these equations of motion when the constraint is applied to ev- 
ery particle, rather than system in general, such as (3). While it is possible that similar equations 
can be derived by Gauss' principle, we believe that the utilization of geometric methods developed 
here leads to important consequences, which are hard to obtain using traditional approaches. 

Second, and the most important contribution of the paper lies in the proof that the kinetic 
theory for our nonholonomic theory can be taken initially to be concentrated on the constraint set 
in phase space, and it will preserve this structure for all times. This observation will allow us to 
derive the nonholonomic kinetic equation (56) which we believe is the main result of this paper. 

Finally, the remainder of the paper is devoted to the analysis of the derived nonholonomic 
kinetic theory. More precisely, in Sections 4 and 5 we find an explicit solution of the full kinetic 
equation, which is inspired by Poiseuille-type fiows. We also show that while the momentum 
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is not conserved, so the traditional fluid approach based on the density, momentum and energy 
conservation laws is possible in our system, we can still derive exact conservation laws that follow 
from the kinetic equations. We conclude by deriving a hydrodynamic model from a cold plasma 
closure of the moment hierarchy. 

The theory we have developed will be directly applicable to all systems having constraints that 
are linear in velocities. This restriction follows from the limitations of the Lagrange-d'Alembert's 
principle. The question immediately arises whether this theory can be generalized to include 
more general constraints, that are nonlinear in velocities and are applied to every particle in the 
ensemble. As far as we are aware, no such theory has been developed, although we beheve that 
the derivation of equation is possible based on the Gauss' principle. The crucial question for our 
paper lies in the possibility of concentrating the density in the phase space on the constraint set, 
or distribution as it is commonly addressed in the nonholonomic mechanics. ^ 

A small digression into the geometry of nonholonomic constrains is warranted here. The main 
difficulty posed by single-particle nonholonomic dynamics is that it lies on a distribution. That is, 
single-particle nonholonomic dynamics takes place on only a certain subset of the position-velocity 
phase space [ ]. In contrast, Vlasov kinetic theory takes place in phase space with independent 
coordinates {x,v), in which one should not confuse x with v. The constraint for the particles may 
be nonholonomic, but this does not imply a relation between x and v. Instead, the nonholonomic 
constraint is imposed in phase space by assuming the probability density is defined on the whole 
phase space, although it is supported on the distribution on which the nonholonomic relation holds 
when we set i; = f . 

In particular, we shall prove that the probability density that is initially concentrated on 
the distribution, will remain concentrated on the distribution for all times. This is the crucial 
point of this paper. The work here shows that this claim holds for arbitrary constraints that 
are linear in velocities. We also believe that it also holds for more general constraints, as was 
demonstrated recently in the context of kinetic theory for the Vicsek model in mathematical 
biology [21, 25, 26]. There, it was shown that the nonlinear nonholonomic constraint |vp = 1 
for every particle preserves the concentration of distribution on the constraint. The open question 
that remains is the following. What is the most general form of nonlinear constraint that has that 
property? We do not know the answer to this question, which is very interesting, but is beyond 
the scope of this paper. 

The system we study here illustrates the challenges that may lie ahead in developing non- 
holonomic statistical mechanics. For example, the system of interacting particles is closed in the 
thermodynamic sense, as there is no energy exchange with the substrate. However, the system is 
not isolated, in the sense that its momentum is not conserved. In fact, neither the momentum of 
each particle, nor the total momentum of the system is conserved. In some particular cases, there 
are additional conservation laws for individual particles (Routh and Jellett integrals) and these 
are useful, as we will discuss below. As has been shown in numerical simulations [20], because of 
the coupling between the rotational and translational degrees of freedom, nonholonomic dynamics 
breaks both the ergodicity hypothesis and equipartition of energy between different degrees of 
freedom, and the definition and treatment of these concepts on the constraint distribution become 

^We hope no confusion arises between the unfortunate colhsion of terms "distribution" as defining the nonholo- 
nomic constrained set, and the "distribution" as defining solutions in generalized function space, e.g., (5-functions. 
These terms are both well-established in the appropriate literatures, and we will need to use both meanings of the 
word in this paper. We hope that the particular meaning of this word is clear from the context. When a clear 
distinction is necessary, the kinetic probability distribution will be called "probability density" . 
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difficult. Correspondingly, familiar concepts such as thermodynamic temperature cannot be defined 
easily. In addition, because of the lack of momentum conservation, care must be taken in deriving 
fluid-like continuum mechanics for such systems. 

2 Rolling molecules and collisionless kinetic theories 

2.1 Euler-Poincare dynamics of a single rolling molecule 

We shall start with a brief derivation of the equations for an individual Chaplygin ball system. This 
is a well-known problem, but recalling the derivation here will allow us to introduce some useful 
notation. This derivation relies on the symmetry reduction principle for nonholonomic systems 
first developed in [27]. For more details on this derivation and some recent results, see [6, 9]. 

Consider an unbalanced rolling ball whose center of mass is positioned at Ix ^ away from 
the geometric center in the (body) coordinate frame in the ball, where I is a length, and x G S*^ 
is a unit vector in the body frame. As the ball moves, it undergoes a rotation 71 G 5*0(3), and a 
translation x G M^. In particular, its geometric center is translated by x, and its center of mass is 
located by the vector ITZx pointing from its geometric center to its center of mass in the (spatial) 
stationary coordinate frame. Let us also consider an external field E acting on the center of mass of 
the ball. This could be an external electric field, gravity or another external potential force acting 
on the particles whose potential energy is E ■ TZx- One writes the Lagrangian of an individual ball 
as the difference between its kinetic and potential energy: 

L(7^, 7^, X, i) = ^ + i ^ ViA) |7^A|2 dA - E ■ 7^x , (4) 

where A is the position vector of a point in the molecule, which in turn possesses mass density V 
and volume B. Here we shall follow the Euler-Poincare symmetry-reduction principle [28]. For 
this purpose, we transform the Lagrangian into the following reduced spatial variables: 

Director, n = TZx G 5*^; 

Angular velocity in the spatial frame, = TZTZ^ G 50 (3); and 
Tensor of inertia in the spatial frame, j = TZJTl^ G sym(3 x 3) , 

where so (3) denotes the Lie algebra of antisymmetric matrices and J is the symmetric inertia 
tensor J = — Jgp(A)(AA^ — dA. The last of these spatial variables (j) is known as the 

'microinertia tensor' in the theory of micropolar media [29, 30]. In fact, micropolar media provide 
an interesting analogy for the systems considered in this paper and many of the methods used here 
transfer easily to the micropolar setting. One difference, however, is that the quantity n = TZx is 
not strictly a director, since n 7^ — n. As we shall see, this parity invariance under Z2 is broken 
for rolling molecules by both the potential E ■ n and the nonholonomic constraint. 

Remark 2.1 (The hat map). In defining the angular velocity variable v, we use the following 
hat-map correspondence between an antisymmetric matrix d G so (3) and a vector G M'^ with 
components Vi, % = 1, 2, 3, 

^ifc = - (^jkiT^i and i/j = - ^ eijfcZ^jfc , i, j, k,l = 1, 2, 3, (5) 
in which etjk is the completely antisymmetric tensor with €123 = 1. 
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The symmetry-reduced Lagrangian corresponding to (4) for an individual ball may now be 
expressed in spatial coordinates as 

^(i^, j,n,x,x) = ^|xp + ^jV u-E n, (6) 

in which the time derivative is denoted with an over-dot, e.g. d'^/dt =: x. The rolling constraint 
for an individual ball reads 

X = n{rn^z + x) = X (rz + n) := X cr(n) , (7) 

where z is the constant spatial unit vector pointing perpendicular to the substrate and r is the 
radius of the ball. In (7), we have also introduced the notation 

cr(n) := rz + n , (8) 

in which the spatial vector cr(n) points from the contact point of the ball to the position of its 
center of mass at a given time. 

Remark 2.2 (Dimensionality of the rolling constraint). The rolling constraint (7), technically 
speaking, defines the relationship between three-dimensional vectors. However, this contraint only 
contains meaningful information about the motion of either the contact point or geometric center 
of the ball, which are related by being offset by given constant vector. Thus, we shall understand 
this constraint as a relationship in TS0{3) x TM^. Similar considerations will apply later to the 
constraint applied to the motion of an assembly of particles. 

Proposition 2.3. The constrained Euler-Poincare variational principle, i.e., Hamilton's principle 
for the symmetry-reduced Lagrangian, 



rt2 

5 I /(r^, J, n, X, x) dt = 

Jti 



yields the following nonholonomic equations 



d 61 61 
at 6u 6h> 



■ 61 
5j 



61 f d 61 61 , , , 

— X n = -——-— X o- n , (9) 
6n \dt6± 6^' ^ ^ ^ ^ 



f^ + m=0, (10) 

^ + nxzy = 0. (11) 

Remark 2.4. The notation ~Xi = eijkAjk for an arbitrary antisymmetric matrix = —A has 
been introduced in (9). Likewise, because of the hat map in (5) we have n x = —vn. 

Proof. The proof of the proposition is standard for this kind of problem, see e.g., [9], and it uses 
the relation <t = ri obtained from the time derivative of (8). □ 

Corollary 2.5. For the reduced Lagrangian (6), the dynamic equation (9) may be written equiva- 
lently as 

^ (jV) + jV X 1/ + ^ [j, i/i/^] — Exn = — (i>xcr)xcr — (j/xn)xcr. (12) 

Then, upon using the properties of the hat map, one may rewrite this equation as 

jO + cr X {cr X 0) =: (j + aa)^ = Exn + jVxi/ + crx (i/x (i/x n)) . (13) 
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Equation (13) differs in form from tfie standard equations for the rolling ball [6, 9]. This is 
because the standard Chaplygin ball equations are written in the body frame, whereas equations 
(13) have been written in the spatial (fixed) frame, in terms of the time-changing moment of inertia 
governed by (10). Although they have been written in the spatial frame instead of the body frame, 
equations (13) are mathematically equivalent to the standard equations for the Chaplygin ball. 
While previous authors have considered the Chaplygin ball equations in the body frame, the 
spatial frame will be preferable for our applications later with ensembles of rolling balls, despite 
the necessity of considering the inertia tensor as a time- dependent variable. 

As mentioned earlier, additional conservation laws due classically to Routh and Jellett exist 
for an individual rolling ball with symmetry. See e.g. [ ] for references and discussions of these 
classical integrals of motion. The Jellet conservation law, in particular, will play an important role 
here in our considerations of multi-particle dynamics. The Jellet integral 

Qj = ■ '^(n) (14) 

is conserved for an individual Chaplygin ball, provided that two of the ball's inertia tensor's 
eigenvalues in the body (let us call them Ji and I2) are equal in value Ji = I2, and the axis of the 
third inertial eigenvector is coUinear with connecting the center of mass and geometric center. 
In that particular case, there are two more integrals of motion: the energy of an individual ball 
and the Chaplygin (Routh) integral. As was shown by Chaplygin [ ], the presence of these three 
integrals makes the rolling ball equations completely integrable. When there are several rolling 
particles of Chaplygin type interacting through a central potential directed at their centers of 
mass (such as the Lennard- Jones potential), the Chaplygin integral is not conserved for either an 
individual ball or the whole system. As one might expect, the total energy of the entire system 
is conserved. However, more remarkably, the Jellet integral is still preserved for every individual 
ball, in spite of their interactions ["20] . 

The next section introduces the main problems that are encountered when trying to build a 
kinetic theory of rolling molecules. As we shall see, these problems emerge fundamentally from 
the absence of an invariant measure in phase space, arising because the single-particle dynamics 
is not Hamiltonian, except in very special cases. 



2.2 Preliminary considerations in kinetic theory 

In order to describe multiparticle dynamics on phase space (typically TM^ ~ M^), kinetic equations 
govern the dynamics of a probability density f{z,t) on phase space. When collisions are neglected, 
kinetic equations are transport equations of the type 

|^ + div,(/X) = 0, (15) 

where the vector field X usually contains nonlocal terms deriving from the collective interactions 
among particles. In the Lagrangian picture, this means that a kinetic equation is given in charac- 
teristic form as 

^(/,(V'(t))d'^V'W) =0 along V'W = X(V'(t)) (16) 

where 'ify{t) is a characteristic curve in phase space. Typically, il^{t) = (x(t), v(t)) G in ordinary 
position-velocity coordinates, so that X(i/') G and k = Q. Alternatively, given the phase-space 
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curve il^{il^Q,t) (with initial label i/^q) and the initial probability density foi'^l^o) d'^V'o) ^^e forward 
probability density function is given by the Lagrange-to-Euler map 

f{z, t) = I M,) 6{z - V'(Vo> t)) d^V'o • (17) 

This expression recovers the well known Klimontovich particle solution for an initial point particle 
/o(i/?o) = ^(V'o ~ ^o) Eulerian point Zq. Here, the notation ft{il>{t)) with time dependence 

in ft indicated by the subscript t is used in the Lagrangian representation, while f{z,t) represents 
the corresponding Eulerian quantity. In mathematical terms, one says that ft is obtained as the 
push-forward of /o by il^{t) and one writes, e.g., 

ft = ^Mk. (18) 

In coUisionless systems (i.e. systems of the type (15)), the motion is reversible. When this 
motion is also Hamiltonian, the flow preserves the Liouville volume on phase space d'^t/jg, so that 
d^il^{t) = d^ij^Q and hence the relation (16) in the form 

ft{m)d'm = fo{i^o)d'i^o 

implies that the probability density evolves as a scalar function. That is, ft{il^{t)) = /o(i/'o) 
Hamiltonian systems. The invariant Liouville measure is the basic ingredient of modern ergodic 
theory. In the thermodynamics of closed Hamiltonian systems, the combination of reversibility 
and the existence of an invariant measure implies conservation of the entropy functional 

S = Jfhgfd'z (19) 

(in units of Boltzmann constant). At this point, one may inquire into the definition and dynamics 
of entropy. As we shall discuss in Remark 4.6 below, the concept of entropy for nonholonomic 
systems produces divergences that need to be treated by mathematical methods that are beyond 
the scope of the present paper. 

The issue of thermodynamic entropy may also be related to the closedness of rolling particle 
systems. For the special case of a body rolling on a substrate, nonholonomic dynamics emerges 
from the constraint force that the substrate exerts on each body. Thus, although this force does 
no work, according to the Lagrange-D'Alembert principle [6], the constraint still represents an 
interaction with the substrate. According to the ordinary thermodynamic definition of a closed 
system, the system would be still be defined as closed, since no energy transfer occurs from or into 
the substrate (at least in the absence of sliding). However, the point is that the conserved energy 
of rolling bodies does not imply Liouville volume conservation and thus ordinary thermodynamic 
considerations do not apply in general to a nonholonomic system. 

Remark 2.6 (Probability density and nonholonomic distributions). As briefly mentioned in the 
previous section, nonholonomic systems are defined on special geometric objects known as distri- 
butions. These are hyperplanes in phase-space to which the dynamics is confined. In the general 
case, distributions possess no differentiable structure and the concept of a probability density may 
not make sense in terms of a differential form of maximal degree. Our strategy in this paper will 
be to define a weak density on the whole phase space, that is supported on the nonholonomic 
distribution. In this way, one can still use the differentiable structure of the entire phase space, 
while particle dynamics is still localized on the distribution. However, any statements about the 
probability distribution will only hold in the weak sense. That is, they will hold when integrated 
against smooth functions on phase space. 
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3 Lagrangian trajectories of rolling molecules 

3.1 Relabeling symmetry of collisionless kinetic theories 

It is clear from (17) that all the information contained in a collisionless kinetic equation is encoded 
in the phase-space transformation i/j(i/7Q,t) taking the Lagrangian label i/^g to its phase-space 
position at time t. Finding this transformation is equivalent to solving the equations of particle 
motion t/? = X(i/') with initial condition t/'(0) = i/'o- Thus, one is normally interested in finding 
the vector field 



which produces the equations of motion. (Here, 'o' denotes composition from the right.) The 
equations of motion for nonholonomic systems may be derived from Hamilton's variational principle 
[9], provided one evaluates on constraint surfaces only after taking variational derivatives. So, for 
non-interacting particles the vector field is well known and it can be derived from constrained Euler- 
Lagrange equations. When particles interact mutually with each other, the Klimontovich method 
[' ] may be applied to derive a collisionless kinetic equation. However, the application of the 
Klimontovich method on nonholonomic distributions may present difficulties, if the differentiable 
structure were to be lost, so that the divergence in (15) would not make sense, see Remark 2.6. 
Even in the case when nonholonomic constraints in lead to an appropriate differential structure 
(e.g., as defined by the implicit function theorem), it can still be difficult to find a convenient 
coordinate system. 

Thus, our strategy in formulating a consistent kinetic theory for nonholonomically constrained 
particles will be first to use Hamilton's variational principle to derive the equations for the La- 
grangian trajectories il^{i(^o, t) on the whole phase space and then to use reversibility to obtain the 
Eulerian vector field X(2) = o jp^^^z). The general formula for the variational principle reads 



which in turn yields constrained Euler-Lagrange equations through the nonholonomic constraint 
in expressing 6i/). The labels i/^q are integrated over the probability density /o(i/'o) d^'^o whose 
support will later be confined to the nonholonomic distribution by using a Dirac delta function. 
This Lagrangian variational approach is widely used in the theory of inviscid fiuid fiows, and its 
application to kinetic equations was first due to Low [23], who successfully showed that Vlasov-type 
equations in plasma physics possess a Lagrangian variational formulation. Forty years later. Low's 
work was revisited in [ ], within the modern mathematical language of Euler-Poincare reduction 
by symmetry. Since the action principle is invariant under relabeling, applying the inverse V?"^ of 
i/' yields the following variational principle in terms of purely Eulerian variables: 



X = tl^ o 



1 




(20) 



= 6 








(21) 
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where / is given by the above mentioned Lagrange-to-Euler map. In the simplest case of holonomic 
particle dynamics, the variations 

5X = dt{6il) o V"^) + o V"^) • V)X - (X ■ V)(5V o '0"^) 
5f = -divif 611^01/^-^) 

yield the equations of motion 

||-+div(/X) = (23) 

Although the equations appear to be coupled, the special nature of the Lagrangian that we shall 
choose will decouple the above system so that the second line acquires its own meaning as a kinetic 
equation. The present paper extends this method to nonholonomic systems, in considering the 
special example of interacting rolling balls. 



3.2 Configuration space and Lagrangian 

The main difficulty in developing the kinetic theory of nonholonomic systems can be explained 
as follows. Since Chaplygin's ball is a nonholonomic system, its dynamics takes place on the 
distribution V C T5'0(3) x TM^ formed by the nonholonomic constraint 

P={(x,v,7^,^;7^)GT(M2x50(3)) : v = ^;7^7^^(rz + Te^)} , (24) 

where TQ denotes the tangent bundle (i.e. the position-velocity phase space) associated to the 
configuration manifold Q, so that TM'^ = M?^ and (x, x) G 750(3) for any trajectory x(t) G 50(3). 
In general, the distribution, (or set), defined by the nonholonomic constraint is not a manifold, 
and it does not possess the familiar tangent bundle structure common in mechanics [6]. Thus, a 
great care must be taken to introduce any familiar concepts borrowed from calculus on manifolds. 
These difficulties persist if we wish to make a kinetic theory of nonholonomic systems. Formally, the 
Lagrangian derivation of a collisionless kinetic theory should involve smooth invertible coordinate 
transformations 1/? (diffeomorphisms) on V, which in turn would determine the Lagrange-to-Euler 
map (17) associated to a probability density / on V. However, restricting the process to respect 
the geometric structure of the nonholonomic distribution V C TM^ x rS'0(3) leads to analytical 
difficulties. We choose to circumvent some of these difficulties by considering dynamics on the full 
phase space TM^ x TS'0(3) for as long as possible, before inserting the constraints using delta 
functions. 

Of course, physically, the dynamics makes no sense outside the distribution V. Thus, in what 
follows, we consider kinetic probability densities that are defined everywhere on TM^ x TS'0(3), 
but are concentrated on the distribution V and vanish outside V. To justify this assumption, we 
will show that kinetic densities concentrated on V remain concentrated on V for all times under 
the evolutionary equations we shall derive. Thus, we define the 4-coordinate Lagrangian mapping 
(diffeomorphism) 

:= (^x,^v,^7e,^.^) e Diff(TM2 x rS0(3)) (25) 

which takes the initial coordinates (xq, vq, TZq, v-jiq) to their values at the time t. Here, Diff (M ) de- 
notes the Lie group of diffeomorphisms of a manifold M. It may be worth emphasizing that, 
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despite our notation, -0 is not simply a vector in 3D; rather it is a Lagrangian map on the 
phase space TM^ x TS'0(3). In terms of this Lagrangian map, the rolling constraint identi- 
fies the following infinite-dimensional nonholonomic distribution in the ambient tangent bundle 
TDiff(TM2 X TSO{?>)): 



V 



|TDiff(TI 



(26) 



In this section, we shall show how the following Lagrangian of the type (20) produces nonholo- 
nomic rolling dynamics: 



/o(xo,7^o,vo,W7^o) 



+ 



+ 



Viz 



dxodvod^;7^od7^o, (27) 



where we have used the relation ip-ji-i = ipn II ' II j denotes the norm given by the trace as 

\\A\\y.= TT{A^JA) . 

At this point the right trivialization map [ ] gives the change of coordinates 

(7^, z/) := (7^, vjin-^) e S0{3) x so(3) . 

where u = v-jiRr^ is an antisymmetric matrix belonging to the Lie algebra so (3) of the rota- 
tion group 5*0(3). Notice that, for ease of notation, we have dropped the hat symbol ^ usually 
accompanying antisymmetric matrices. See Remark 5. We then have Diff(T]R^ x TS'0(3)) = 
Diff(S'0(3) X so(3) X TM^), so that the new Lagrangian reads 



/o(xo,vo,7^o,z/o) 



- 2E ■ i^nX 



+ 



dxodvodz/odT^o 



(28) 



To summarize, we start with the tangent space TDiS{TQ) of diffeomorphisms of a tangent bun- 
dle TQ = T{S0{3) X M^) (not a distribution). Eventually, these act on a probability density 
/o € Den{TQ) that is also defined on the same tangent bundle TQ. Then, enforcing the nonholo- 
nomic constraint ?/'x = 4'n {r'lp-^ z + x) on the diffeomorphisms yields a dynamics that lies on the 
constraint distribution V C TDiff(TQ). So far, this is the standard symmetry reduction approach 
to nonholonomic systems [ ], although now it is being applied in infinite dimensions. 

The key point of this paper comes from the fact that the dynamics of the individual microscopic 
particles take place on the constraint distribution V C TQ. However, since V is not even a 
manifold, we cannot consider such an object because there is no sense in which a probability 
density is defined on the constraint distribution - at least, we are not aware of such work and in 
any case, the consideration of functions on V leads to severe analytical difficulties. Instead, we 
avoid these difficulties by introducing a natural construction that is the key step in this paper. 
Namely, we require that the probability density /o G Den{TQ) defined on the whole tangent bundle 
TQ be supported on the distribution V C TQ. This is done by taking the following singular ansatz 
for the probability density 



/o(xo, vo, Uq, 7^o) = 0o(xo, 1^0, 7^o)5(vo -i^ox c^(7^o)), 



(29) 
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then showing that the dynamics of the resulting kinetic theory preserves this class of solutions. 
As we shall see, this singular ansatz for the probability density leads to Euler-Lagrange equations 
involving the whole tangent bundle TQ, which are supported on the nonholonomic distribution 
V C TQ on which the particles are constrained to move. This is a natural picture, since the 
dynamics does not "see" what is outside the distribution V. In this way, we avoid the difficulties of 
dealing with the densities on V only and, as we show below, derivation of kinetic theory is possible. 
We believe that a similar approach can be generalized for an arbitrary system of nonholonomically 
constrained particles, without many substantial difficulties. Notice that the probability density 
(po = J fo d^Vo is not a density on the distribution V] rather, this is more simply a probability 
density on x TS0{3) ~ ^ S0{3) x so(3). 

3.3 Hamilton's principle and Euler-Lagrange equations 

The Euler Lagrange equations arise from Hamilton's principle 

6 r Lf,{^,^)dt^O. 
Jtl 

The variations are subject to the constraint 

(^V'x = ^ipniPn (^i^n) (30) 

Then, upon denoting the Liouville volume by dwo — dxodvodi/odT^-o, one computes the following 
relation 

54 + ^ • Si/j^]dwodt = - / ( f ^ - ^ ) • Stp^dwodt 



ti 



d SLf^ 




dt 5tp^ 




dSLf, 




dt 5ip^ 





(Si'n^Pnf^ii'n)) dwo dt 



-i:k 



where the superscript ( • denotes antisymmetric part. Next, Hamilton's principle yields the 
Euler-Lagrange equations 



dtSi^n S^nJ \dt dijj^ ^V^x / / 

0, (31) 



d 




dt 6ip, 




dSLf, 




dt Sij^ 





= 0. 



The last two Euler-Lagrange equations lead to the relationships 

/o (V'x - V'v) = and fo (j^n - '^u'^nj = , 
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where we keep in mind the singular expression of /o given in (29), so we do not divide out by 
/o- Then, upon introducing cri^pTi) = rz + ipnXi the nonholonomic constraint ipy^ = ip-jitpj^crlipTi) 
yields 

/o(^v-^. xo-(^7e)) = 0. (32) 
Notice that the presence of a delta function in the density /o forbids dividing by /q. 

Theorem 3.1 (Nonholonomic Euler-Lagrange equations). Subject to the constraint (30), Hamil- 
ton's variational principle associated to the Lagrangian functional (28) yields 

fo (^x - X cr(V^7e) j = and fo (Jj-r - ipuipnj = , 

Also, upon introducing the director nlipji) = ipnX ^'^^ the microinertia tensor jlip-ji) = ip-jiJip'^, 
the resulting nonholonomic Euler-Lagrange equation reads 

fo[j{^n) ipu - Ji^n)^u X Viv) = /o(e X niipn) - criipn) x (T{ipn) x i^u - (T{ipn) x cT{iljn) x 

(33) 

Proof. The first part of the theorem has already been proven above. We only need to prove the 
second part. After computing 



dtp' 



'7^ 



^ = /o((Ex^^D-^ + ^b,^.V^.])V^7e, 



6Lf^ 



0, 



the remaining equation gives the dynamics of the local angular momentum 

n:=^V^^G50*(3). 
Indeed, pairing the first equation in (31) against 6tpTi and differentiating by parts yields 



^^^^ ,h \ +( ,lr'\ ((^ ^^f' - ^ 



V v 

micropolar terms nonholonomic terms 



Therefore, the desired equation of motion is 



A 



fo ^(J(^7^)^.)'' = - /o (J(^7^)^. ^P.)^ + fo (ex^^^ + ^ [j, ^,^.] ) - /o o-^)-" 
= /o(En^(^^))^-/o(V'v^^)^. 
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Since dj('?/'7e)/dt = [^puij], we then find (33), upon using the inverse of the hat map and its 
properties (see e.g. [■)■ ])• ■ 

Notice that, upon dividing (33) by the 0o(xo, i^o, "^o) in (29) and integrating over vq, we obtain 



{ji'u — jil'v xipu — 'E'Xn — (Txcrxijj,^ — (TX&xipL 



= 0. 

vo=i^ox<t(''^o) 



Consequently, Lagrangian dynamics is supported on the distribution (24), consistently with the 
single molecule dynamics. Moreover, this calculation shows that the ansatz (29) remains consistent 
for all times. 



4 The Euler-Poincare approach to kinetic theory 

As an alternative to the Euler-Lagrange formulation, let us derive these equations of motion using 
the Euler-Poincare theory, which explicitly utilizes the symmetry reduction in the Lagrangian. 
This reduction arises from the relabeling symmetry introduced in equation (21). 



4.1 Nonholonomic Euler-Poincare equations of motion 

This section makes use of the following invariance relation (see (21)): 

Lf,{^P,^P)=L^,f,{^Po^p-') =:1{XJ) 

where X = ip o tp^^ is the vector field transporting the probability density / = i/'*/o) where ^/'^./o 
is the push- forward given by the Lagrange-to-Euler map (17) 

/ = y foiwo) 5(x - ^x(wo)) <5(v - V^v(wo)) 5in - ^7^(wo)) S^u - i!u{wo))dwo 

= j t/'o(ao) (^(x - V'x(ao)) ^(^ - ^7^(ao)) 5{v - ^^.(ao)) 5(v - ^v(ao)) doo . (34) 

Here we have introduced the notation ^'('^o) '■= i^{'^o)\wo=uo^tT{'Ra) 

Wo := (xo,Vo,7^o,^'o) , := (xq, 7^o, i^o)- 

Taking the time derivative of the Lagrange-to-Euler map (and pairing it with a test function) 
produces the evolution equation for the probability density 

^ + V-{fX) = (35) 

where X = X(x, v,i/,7^) = o 

Upon using the Lie derivatve notation £x (see [2<S]) and by introducing rj = Sip o the 
variation 6X = rj + [X, rj] yields 

/•t2 rt'i 

6 l{X,f)dt= / 
Jti Jti 



'dtJx'^^'Jx^^^J]-'^ 



dt. 
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where 51 /5X is a differential one-form density on x5o(3) x 50(3) (see [ ] for tlie explicit forms 
of tfie Lie derivative) and (■, ■) denotes the pairing between vector fields and one-form densities on 
the same space. Then, we make use of the constraint 

'7x = ?77^7^^<T(7^) , 

and the Euler-Poincare equations read 

d 51 ^51 f d 51 ^51 ^„5l' 



ordinary terms nonholonomic terms 

d SI 51 „«\ „ i d St „ 51. ,„5t . 



Following [20] , in order to account for collective interactions among molecules, such as dipole and 
Lennard- Jones, we also introduce an appropriate interaction potential U = U{'x — x.',TZ'TZ'^), which 
generates a term 

U* y"/dvdi/ = j t/(x-x',7^'7^^) y"/(x',v,^y,7^')dvd^/dx'd7^' 

in the Lagrangian. Then, the Euler-Poincare Lagrangian reads 



+ \u^-^f + ||u7^7^^ - ^^||^^dxdvdI/d7^. (38) 

The other two equations of motion yield 

/(n7^-z/7^) = 0, / (n, - v) = , (39) 

so that finally 

fx = f {■v,un,a^,a^) . 

Notice that upon multiplying the nonholonomic constraint Mx = u-jiTiT- criJZ) by /, the equations 
above would yield 

/(v-j/ X ^T(7^)) = 0. (40) 

If / were a smooth nonvanishing function, this relation would be contradictory, because it would 
imply a relation between independent coordinates. However, the singular expression for /o in 
equation (29) will be seen to avoid this difficulty, as a result of the following Lemma, which will 
be proven in two different ways. 

Lemma 4.1 (Singular probability density function). If f^ is given by (29), the definition f = il'*fo 
and the relation (32) imply 

/(x, V, u, 7^, t) = 0(x, u, n, t) 5{y -V X (rijl)) (41) 

where 

0(x,^/,7^) = (^/i*0o)(x,^/,7^) = / 0o(ao) "^(x - V^x(ao)) ^("^ - ^7e(ao)) ^('^ - V^i.(ao)) c^cto • (42) 
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Proof. Recall the relation (32) fo[tpv — x (^) = 0- Integrating this identity over vq yields 

V'v(ao)) = ijAao)) X cr(V^7^(ao)) . 
Then, the Lagrange-to-Euler map (34) becomes 

/ =y0o(ao)'^(x- V^x(ao))^('^-^7^(ao))^(^^-^I^(ao))^^(v- V^^(ao) x cr(^7^(ao))) dao 
and thus taking the pairing with an arbitrary test function /i(x, v, i/, TZ), we have 
J /i(x, V, u, TZ) /(x, V, I/, TZ, t) dxdvdz/dT^ = 

= j h{'i(j^{ao),4i^{ao)) x c^(v5^7^(ao)), V^i^(ao), V^7^(ao)) 0o(ao) dao 
= jdxdudn jdv /i(x, v, u, TZ) 0(x, u, TZ) 5{\ - u x CT{n)) 
= jd\ j dxd^/d7^ /i(x, v, i/, TZ) 0(x, z/, 7^) (5(v - x cr(7e)) 
/;,(x, V, i/, 7?.) 0(x, z/, 7^) (5(v — x (t{TZ)) dxdvdz/d7^, 
where the intermediate steps follow by direct verification. Consequently, 

jh{x, V, u, 7^) (^/(x, V, z/, 7^) - u, TZ) 5{v - u x c^(7^)) j dxdvdz/d7^ = 
and the statement of the Lemma follows, since h is arbitrary. ■ 

Remark 4.2. Thus, the potentially troublesome equation (40) presents no difficulty, as it is sat- 
isfied automatically because the kinetic probability distribution f takes the form (41)- Namely, 
because f takes the form (41), equation (40) follows from the relation for delta functions that one 
interprets xS{x) = for argument x. In the next section, we shall show an alternative proof of this 
Lemma that utilizes the dynamics of the kinetic equation. The interpretation of equation (39) will 
be discussed further in Remark 4-5. 

At this point, one may proceed by deriving the Euler-Poincare equations from the following theo- 
rem, whose proof may be found in Appendix A. 

Theorem 4.3. Upon using the reduced Lagrangian (38), the nonholonomic Euler-Poincare equa- 
tions (37) produce 

f {un - vTZ) = , / (ux - z^ X cr) = . 

Also, define the director n(7^) = IZx and the microinertia tensor j(7Z) = IZJIZ^^ . Then, upon 
recalling the notation := tijk Ajk for an arbitrary antisymmetric matrix A, the Euler-Poincare 
equation (36) yields 

f {j0.u + (T X (T X at,) 

= f ^ju xiy + Exn - (^dnU *jf dvdi/^ 7^^ + c^x(^yx^/xn) + c^x d^U * / j . (43) 
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Therefore, the equation for the i/-component of the Euler-Poincare vector field X reads 
fcbu = f {j + (^(^y^ X i/ + E X n — ^TiU * J f dvdi^-^ TZ^ + crx(i/xi/xn) + crx d^U * f 

where aa = crcr^ — cr^I is a traceless symmetric matrix that is produced by the nonholonomic 
constraint and modifies the microinertia tensor j. 

4.2 The probability density function 

Upon recalhng the Lagrangian dynamics arising from Euler-Lagrange equations, the Lagrange-to- 
Euler map / = t/'^k/o in (17) can be apphed to give the exphcit expression of the probabihty density, 
as was done in the previous section. Then, the time derivative of the Lagrange-to-Euler map (ap- 
propriately paired with a test function) yields the kinetic equation (35) with X = X(mx, uh, a^, ay,). 
This kinetic equation becomes 

% + ^^- (/^x) + Vv ■ (/av) + V, ■ (/a,) + Vn ■ ifun) = (45) 
where the relations 

fu-iz = fvTl and fu^ = x c^(7^) 

were found in the previous section, along with the expression (44) for fai^. At this point, upon 
recalling the relations (18) and (32), we observe that /ov can be computed from the Euler-Lagrange 
equations as follows 

/av = /^v o^-^ =f {^^,(T{il)n) + iJuCriipn)) tp'^ 

= /V'.rV(7^) + fu (^7e^"i ■ dn)(T{n) 
= fa^cr(n) + fv [un ■ 57^)c^(7^) 

= fa„ X criJZ) + X ly X n. (46) 
Consequently, the above expression for av allows to write the kinetic equation (45) as 



This result allows us to demonstrate an alternative proof of Lemma 4.1, by writing the evolution 
equation in the weak form. 

Lemma 4.4 (Solutions on the distribution). Suppose the initial condition for the probability func- 
tion f are of the following form 

fit = 0, X, jy, 7^, v) = 0o(x, 1^, 7^)(5(v - x c^(7^)) , (48) 

Then, at arbitrary t > 0, 

f{t, X, u, 7^, v) = (f){t, X, u, 7^)5(v -ly X c^(7^)) , (49) 
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i.e., the concentrated solution preserves its form under evolution. Moreover, evolution of (p is given 
by the equation 



^ + ^/xc^•|^-Tr(7^^z7 
at ax 



d 



(50) 



Proof. Consider, for now, / = (f){^,i',TZ,t)g{\ — x criJZ)), where g is an arbitrary generalized 
function. Then, for an arbitrary test function C(v), we have 



Vdt 
+ 



+ J/ X cr ■ ^ - Trin' u 
ox. \ 



+ 



d 



\ dg dg dg 

on Ou avJ 



8nr eV^^^A^^^^^). 

,C(v) 



0. 



(51) 
(52) 



The pairing, denoted by <, >v, is assumed to be the usual pairing in v coordinate only. Then, 
since only g depends on v, and g = (5(v — i/ x cr), we can rewrite (51,52) as 



90 



— a,, X (T — u X u X n 



d 



on J du 



1/ X cr ) = . 



(53) 
(54) 



From the definition of cr = i; + TZx = z + n, with x being a constant vector denoting the distance 
from the geometric center to the center of mass in the body frame, we get 



d 



a 



dn 



= [Uxf = n^. 



The expression in square brackets in (54) was obtained using the following useful identities for the 
derivatives of 

Tr(^^7^^) , r/(v)|) = (^Tr(z>^i)n) , C(v)^ = ( - ^ x ^ x n) ■ |^(^ x a) . 



Also, 



and 



a.-^,C(v) 



av ■ 1^ , C(v) 



«v ■ ^{i' X cr) 



Therefore, the equation in square brackets in (54) is equal to 



a,, xcr — lyxi/xn 



^^,uxct)=0, 



according to (46). 

From the coefficient of ({u x cr), we see that the equation (50) is valid. The coefficient of VvC) 
evaluated at the point u x cr, vanishes identically due to (46). Thus, the solution (49) remains 
concentrated on the constraint distribution v = x cr(7^) for all times. ■ 
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Remark 4.5. A similar calculation performed in the strong sense, without multiplication by a test 
function C(v), proves that the structure of any solution of the form 



f = ly, Ti, t)g{-v -u X (t) 



(55) 



is preserved under the evolution for an arbitrary smooth function of g. However, we cannot assign 
any physical meaning to this interesting mathematical fact, as the solutions of the type (55) do 
not concentrate the probability function f onto the constraint manifold. Alternatively, if one were 
to take g{$,) ^ in (55) above, then relation (40) would make no mathematical sense and the 
present construction would not be consistent. 

Remark 4.6 (Entropy). Equation (47) produces the dynamics for an entropy 



S 



f\ogfd?xd?ud^nd?v 



where d?TZ is the Haar measure on S'0(3). However, given the singular form of f in (49), this 
expression poses severe problems in the definition of Gibbs entropy S and in its dynamics. Indeed, 
one can use (35) and (49) to write 



dS 
dt 



- fdivX d^x dV d^n d\ = - 



;divX)|,^,,^(^)d=^xdVd37^, 



Unfortunately, the quantity {div X)\^^^^^^^^ is undetermined; since all that is known from (49) 
is the value X\v=uxcr{n)- ^^'^^ point, one may be tempted to require (divX)!^^^^^^^-^^ = 
for physical reasons. However, the exact mathematical (and even physical) interpretation of this 
condition is not clear to us, thus we we shall not discuss this possibility further here. Notice that 
this indeterminacy problem is not related to the compressibility of the phase space flow and thus it 
cannot be overcome by the insertion of a metric tensor, as in equation (2). 

A more convenient form of writing equation (50) may be obtained by introducing the density 
variable 

V5(x, V, j, n, t) := 0(x, u, 7^, t) 
where j = TZJTZ^ and n = TZx- Taking the differential of the above definition yields 



dn 

so that the kinetic equation reads as 



dip r 
on 



A 



dip 
dj 



,J 




The above equation emerges from the moment dynamics for equation (47) and incorporates its 
physical content. Strictly speaking, the moment equation (56) cannot be interpreted as a kinetic 
equation on its own, because the space on which it is defined is not an appropriate phase space 
in position- velocity coordinates (there is no velocity associated to the spatial coordinate x). In 
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particular, a relevant difficulty emerges when one tries to redefine Gibbs entropy in terms of the 
variable alone. Namely, conservation of j \p log d^x d'^z/ d^j d^n does not hold and one is forced 
to adopt ad hoc methods such as the metric tensor approach [15] outlined in the Introduction; see 
equation (2). However, the explicit implementation of the metric tensor approach on a space that 
is not a phase space is beyond the scope of this paper. We just note here that in contrast to these 
earlier works we failed to find an explicit expression for the metric tensor for our case. While it 
may be possible that such an expression could be found, we have left this interesting exploration 
to future studies. 

Here, we assume f/(x — x', IZ'IZ^) = W(x — x', n ■ n'), so that 
a„(x. V. n.,) = (, + 55)- X . + E X „ - „ X 3M d.d, 

+cr x(i/xi/xn) + crx * j di^dj^ (57) 

where we write, for brevity in notation, 

U *j^A.i, = /m(x - x',n . n') ^(x. n.,) d.d, dx'dn' . 

Remark 4.7 (Individual particle solutions). One may verify that equation (47) admits the (Klimon- 
tovich) single-particle solutions of the form 

ip = 6{^- X{t))6{iy - V{t))S{n - N{t))S{j - Jit)) (58) 

where X{t), V(t), N(t), J7'(t) satisfy the single particle solutions for the individual ball (9, 11, 
12), and the rolling constraint 

X{t) = V X o-(N). 
4.3 Conservation of the Jellett quantity 

Let us now move away from considering an individual particle, and move towards the assembly of 
particle given by the Vlasov equation. Since the motion of individual particle has particular con- 
servation laws (under proper assumptions corresponding to Chaplygin integrals), it is interesting 
to see the modification of these conservation laws in the kinetic framework. 

In what follows, we shall concentrate on the microscopic Jellett quantity 

q, = [ju) ■ <T(n) . (59) 

As discussed in Section 2.1, the Jellet conservation law plays an important role in our considerations 
of a kinetic theory for multi-particle dynamics. For a microscopic particle, we write the following 
relation 

^ = j^{3u, o-(n)) = Qu + ju , o-(n)) + {ju, ^ri) = (60) 

due to the conservation of the Jellett integral for an individual ball. In terms of the kinetic 
approach, equation (60) is reformulated as the identity 

( - b, + 3^u) ■ o-(n) - ju ■ (— ■ J/ X n) = , (61) 
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in which ai, replaces i> which is given by (57) and we have also used equations (10,11) to substitute 
for j and n, respectively. It is more illuminating to write this expression as 



dn -X" + ^-^ + Tr^b-,P]^j=0, (62) 

as it is true for any kinetic generalization of a conserved quantity qj. This quantity is essentially 
a generalization of a microscopic time derivative. We define the local Jellett density Qj(x, t) as 

Qj{x,t) = j qj{j,iy,n)(f{x,t,iy,n,j)diydndj (63) 

and compute its rate of change as follows: 

dQj f 9v9(x,t,jy,n,j) , 

-df = J dt ^"^"^^ 

f / dip dip ^ f, d(p\ d . X \ , 1 1 

= -y X ■ ^ + ^ X ^ -Tr(^[j,z.] -J + - . (v.a.) jd^dndj 

J [qjipiy X cr^dudndj + J qjip {diVn(^' x n) + divj[j, i>]} 



d_ 
(9x 



+ /" ^ ■ X n + 1^ ■ + Tr f [j, z>] 
J idn dw \ 



dj 



di^dndj 



Since the term in square brackets on the right-hand side vanishes due to (61), and the divergencies 
inside the curly bracket terms vanish because of the antisymmetry conditions, we conclude that 



dQj _ _ d 
dt dx. 



J (^qjipu X cr)dvdndj . (64) 



Remark 4.8. Note that an exact conservation law for the Jellett quantity no longer exists in the 
kinetic framework. This is because in kinetic theory each point in space contains many particles. 
While the Jellett integral is conserved for each individual particle, the most general conclusion one 
can reach for an assembly of particles is the continuity equation (64). 

Remark 4.9 (Jellett integral for individual particle solutions). One can verify that for the indi- 
vidual particle solutions given by (58), the Jellett integral is conserved exactly, i.e. 

dQAX{t)A) ^ dX 

^ ' = along — := V x <t(N) . (65) 

This is the conservation law of a 3-form density (volume) in the Lagrangian coordinates. Here the 
notation follows directly from (58). 

Let us extend this result for a general conservation law that is satisfied by an individual particle. 
As we have mentioned before, if the particles do not interact, and satisfy Chaplygin's conditions 
formulated above, then each particle has three conservation laws: energy q^, Jellett qj and Chap- 
lygin (or Routh) integral g^- In general, if there is a conservation law for an individual particle, it 
generalizes to the continuum conservation law as follows. 



Holm, Putkaradze and Tronci 



Collisionless kinetic theory of rolling molecules 



22 



Theorem 4.10 (Conservation laws for nonholonomic kinetics). Suppose q{j,iy,n) is a conserved 
quantity for the motion of individual particle, i.e. dq/dt = when v satisfies (11,13). ^ Define 



Q(x,t) = y q{i,v,\i)(^{t,yi,v,rv,i)dudiLvdi . (66) 
Then, Q(x, t) satisfies the continuity equation 

j {qipu X (T{n))dudndj . (67) 
Moreover, Q{'K,t) is conserved exactly on the individual particle solutions (58). 



dQ _ d 
dt (9x 



Proof The proof follows the proof for Jellett conservation law. 

In order to compute the corresponding conservation laws for the energy and Chaplygin ball, 
one defines the corresponding energy density 

Qe(x,t) = j qe{iy,n,j)(p{t,x,u,n,j)dudndj 

and Chaplygin density 



Qc(x, ^) = j Iciv^ n, i)ip{t, X, u, n, j)dudndi . 



This leads to the following Corollary. 

Corollary 4.11. The energy density Qe and Chaplygin integral Qc satisfy the following conserva- 
tion laws 

dQe d 



dt dx 
dQ^ __d 
dt dx 



J (^qeLfu X (T(n))diydndj , (68) 
J (^qc^piy X (T{n))dudndj . (69) 



Note that for the particles possessing arbitrary radial interactions through their centers of mass, 
the Chaplygin integral and energy are not conserved. However, the Jellett integral is still conserved 
for every particle. Consequently, we have paid particular attention to the Jellett integral. 

Remark 4.12. Note that the momentum is not conserved for the system we consider. Neither an 
individual ball, nor the system as a whole conserves momentum, due to the rolling constraint. We 
shall thus avoid deriving equations for the momentum for the fluid of rolling particle, since they 
do not have the physical justification invoked in the motion of "regular" fluid. 



5 Cold fluid closure 

5.1 Cold fluid equations of motion 

Although the conservation laws derived in the previous section may be useful, they are not sufficient 
to close the system and reduce the motion to observables, i.e., to write a fluid-like equations in x 



^It is important that the quantity q does not depend on x. 
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and t. In order to obtain such a closure, we use the traditional moment method for cold plasma. 
In particular, we compute the equations for the moments 

V9 dn dj di/ , 1 ^ dn dj di/ , / n dn dj di/ , j cpdn dj du . 



These equations are found to be 



dt 
d_ 
dt 
d_ 
dt 
d_ 
dt 



J (f dndj du + — ■ J{iy X cr(n)) ip dn dj dz/ = , 
J (fdn dj du + — ■ J{i^ X cr{n))u ip dn dj du — J ai,(pdn dj dz/ = , 
Jmpdn dj dz/ + ■ J{v x cr(n))n ip dndj du — J ip i/ x ndndj du = . 
J j (pdn dj dz/ + ■ J(iy X <T(n))j (p dn dj du — J ip [V, j] dn dj dz/ = . 



The last terms in the third and fourth moment equations arise from integration by parts. At this 
point, one may close the system by invoking the cold-fluid ansatz, 

(^(x, n, J, t) = p(x, t) 6{iy - a;(x, t)) 6{n - n(x, t)) 6{j - J{^, t)) . (70) 

The cold fluid equations follow from this ansatz, dividing the last three moment equations by p. 
The dynamics of the density p is then computed directly, leading to: 

^ + V-{pu:x(T{n)) = 0, (71) 
d(jj 

— + (uj X cr(n) • V)u: = a , (72) 

— + (cj X cr(n) ■ V)n = uj x n , (73) 

^ + (u:xa{n)-V)J=[Q,J] . (74) 

The angular acceleration a on the right-hand side of the moment equation (72) for the evolution 
of spatial angular velocity ct;(x, t) in the cold fluid ansatz (70) is given explicitly by 

a(x, t) = [j + a{n)a{n)) ^ (^Juj x uj — n x dn j p(x') U2{n(x.) ■ n(x')) dx' 

+ E X n + cr(n) x {u: x u; x n + VWi * p) ] ■ (75) 



Here the * symbol denotes convolution and we have chosen a potential Z//(x — x',n ■ n') of the 
additive form 

W(x - x, n ■ n) = Wi(x - X ) + ^2(11 ■ n ) . (76) 

The mass density p affects the angular acceleration in equation (75) only through the interaction 
potentials. Finally, the angular velocity evolution equation (72) may be assembled as 

{J + a{n)a{n)) + (cj x o-(n) ■ V)cj^ = 

Ju) X UJ — n X dn j p(x') Z//2(^(x) ■ 'n.(x')) dx' + E x n + cr(n) x [uj x uj x n + VUi * p) . (77) 
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Equations (71, 73, 74) and (77) provide a closed system. In these equations, the nonholonomic 
constraint from kinetic theory enters equation (72) through the relation for acceleration (75). 

The cold plasma motion equation (77) is formulated in terms of observable variables. One may 
now seek the solutions of the cold plasma closure equations. Of particular interest would be the 
flows in confined geometries, such as straight channel channel. This direction of research will be 
explored in our further studies, as the main focus of this paper is the kinetic theory. Here, we 
present an interesting solution of the full kinetic equations, that is relevant for the cold plasma-like 
flows of the Poiseulle type. 



5.2 Exact solution of kinetic equation (56) in the cold fluid class 

We shall show how to write exact solutions of (56) and (57) in a geometrical setting similar to the 
Poiseulle flow. Namely, let us consider a statistical ensemble of rolling balls, whose direction of 
motion is along the Xi-axis, and whose solution is independent of Xi, but depends on X2. Let us first 
look at the ensemble from the microscopic point of view. We need to enforce that a microscopic 
particle rolls along the xi-axis with a constant speed. This is possible to achieve in Chaplygin case, 
when each particle is symmetric, so ii = 12 (two components of the microinertia tensor are equal), 
and the third axis of inertia is coUinear with the director from the center of mass to the geometric 
center, i.e., x II ^3. If these assumptions about the properties of the microscopic particles are 
satisfied, the rotation about es-axis leaves the tensor of inertia invariant, so j = i. This result is 
afforded by the following 

Lemma 5.1 (Rotation about the inertia axis of symmetry). Suppose TZ is a rotation about the 63 
axis of inertia by the angle a, and the body-frame microscopic tensor of inertia is i = diag{\i, \2, is). 
Then, 



j = 7^i7^"^ = diag(ii, ia, is) + (i2 - ii) 



sin^ a — cos a sin a 
cos a sin a sin^ a 



(78) 





Proof. Proof of this Lemma is obtained by direct computation. □ 

Thus, in the Chaplygin case of a symmetric unbalanced ball, ii = 12 and i = j, so the tensor of 
inertia is the same in the body and spatial frame for this particular motion of rotation about the 63 
axis. In addition, for such motions the position of the center of mass of this particle does not change 
in time, so cr =const, and ^'=const, since the ball will move indefinitely with a constant speed. 
Moreover, since u || 63, it is easy to see that jv || u, so jiy x 1/ = 0. In addition, since the rotation 
is about the x axis, we have n = TZx = X- the absence of external forces and interaction 
potential, all microscopic particles will move independently. Such a microscopic solution sets up a 
reasonable ansatz for the full solution of the kinetic equation. 

Let us now turn our attention to the kinetic equation (56). We need to enforce the microscopic 
motion preserving the Poiseulle flow geometry, so we assume 

(p{t,x,u,n,j) := ^po{x2)5{u - iyo{x2))5{n - no{x2))S{iy - iyoix2))5{j - i) . (79) 

Here, i is a given constant matrix, having the physical meaning of the microscopic inertia matrix 
as described above. Following the microscopic picture, we assume the functions cro(a;2) and iyo{x2) 
to be both perpendicular to the Xi-axis, and \\ n, so ctq x i/q \\ xi-axis. Then, i/q x no = 0, 
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Figure 1: Microscopic realization of Poiseulle flow. Each individual ball rolls in a straight line 
along Xi, rotating about its axis of symmetry 63. The dynamical quantities are allowed to change 
with X2, but since all the balls roll along parallel straight lines without collisions, the solution is 
stationary. 



and since the rotation is only about the third principal axis of the tensor of inertia, one finds that 
i = i and is independent of j. 

The setup of the problem in microscopic setting, is illustrated in Fig 5.2. This leads to the 
simple equation 



d_ 



^0) 



A possible way to satisfy (80) is to posit a^, = 0. Looking at equation (57), we see that Juq x i/q = 0, 
and = iff 

n X (E + dJJ * (p) = (T X dJA * . (81) 

In the case E = (no external forcing) and f/ = (no interaction potential), equation (57) is 
trivially satisfied. 

Another approach for solving equation (80) is to consider U = 0, but E = Eq 7^ 0. The simplest 
condition for solutions to exist is then E || n, i.e., the external field acts along the rotation axis, 
which is also the symmetry axis of the ball. In this case, there is no acceleration of the microscopic 
particle and the equation is trivially satisfied. That case, however, is physically trivial and therefore 
of no interest to us. In a more general case, E need not be parallel to n, even though a^, 7^ 0, since 
the acceleration is independent of u. An even more general solution may be obtained for a central 
potential, i.e. 

W(x-x',n-n') =W(|x-x'|) ■.= U{q). 

In this case, the convolutions in (57) give dJJ * V9 = 0, and we are left with the following integral 
equation for the profile: 



no(x2) X E = (t{x2) X 



x-x' dU 
|x — x'l Oq 



V9o(x2)dxidx2dn' . 



^2) 



One should also remember that cr(x2) = rz + no(x2). It is perhaps easier to solve the equation 
(82) backwards. That is, given a function ip{x2), find E such that the equation is satisfied. 
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Choose an arbitrary function '^o{x2) that satisfies proper smoothness and integration properties. 
For a given W, define 

q(x) = / — v?o(a;2)da;'ida;2dn'. 
J |x — x'l af) 

If we define b = E — q and c = i?z x q, then the equation (82) may be rewritten as 

no(a;2) x b(2;2) = c(x2) , (83) 
with n(a;2) being unknown. The solution to this equation exists if 

b ■ c = (E - q) ■ (z X q) = , 

or simply that 

E ■ (z X q) = , (84) 

which can be satisfied if E is chosen to be perpendicular to the plane of rolling, or E || z. Inci- 
dentally, E II z is the most interesting physical case, describing the case of strong gravitational 
attraction of rolling particles to the substrate on which they roll. 

Given that (84) is satisfied, in order to find the solution for 11(^2) we take the cross-product of 
b(a;2) with equation (83), in order to obtain 

(|bpld-bb^) ■no(x2) = b X c. (85) 

If, for a given b(x2), a homogeneous solution of (85) is given by n/i(x2), and the inhomogeneous 
solution to (85) is ni(x2), then the general solution to (85) is simply 

no(x2) = C{x2)nh{x2) + Ta-i{x2) , 

where C{x2) is an arbitrary scalar function of X2- 

The range of validity for the solution found in this section is narrower than that of the cold fiuid 
closure. However, it affords substantial flexibility in the choice of velocity proflle. This solution 
also illustrates that familiar concepts from the inviscid two-dimensional channel flow, such as the 
choice of an arbitrary velocity proflle, seem to carry over to the nonholonomically constrained 
fluid. In future work, we will address the stability of the derived solutions, which has not been 
addressed here. 

Remark 5.2 (Connection to the cold fluid closure). It is interesting to connect the solution derived 
here to the cold fluid closure obtained earlier in (71-73) and (77). In the assumptions used here, 
we see that both uj and cr lie in the plane perpendicular to the X\-axis, and all variables depend on 
the coordinate X2 only. Thus, x cr ■ V = for all variables. On the other hand, because of the 
choice of the axis of rotation and the symmetry of the microscopic moment of inertia, JT"] = 0. 
Thus, the equations (71-73) are satisfied identically provided that a = 0. Equation (82) is precisely 
that condition with the additional simplification of U2 = 0, so the potential interaction depends 
only on the Euclidian distance between the microscopic particles. 



6 Conclusions 



This paper derives the kinetic theory of an ensemble of interacting particles that are each subject 
to the rolling constraint. The main difference of the present work with the previous literature in 
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the area is that we consider the constraints apphed to every particle in the ensemble, rather than 
on system in general. We have shown how to derive the equations of motion using the geometric 
Euler-Poincare principle, with constraints treated by the Lagrange-d'Alembert principle. The 
limitation to constraints that are linear in the velocities, imposed by the Lagrange-d'Alembert 
principle, may be overcome by using the more general Gauss' principle of least constraint. 

We have not pursued this line of investigation further here, because the Gauss method becomes 
rather cumbersome when the constraint is applied to each particle. In fact, we are not aware for 
any physically meaningful, nonlinear non-holonomic constraints that have been applied to every 
particle, and not to the system as a whole. Nevertheless, it is an interesting topic from the 
mathematical point of view that will be addressed elsewhere. 

The present paper showed that a probability density that is initially concentrated on a con- 
straint distribution that is linear in velocities will remain concentrated on that distribution for all 
times. We are hopeful that this may also be true for constraints that are nonlinear in velocities. 
Recent work on the kinetics of the Vicsek model [24, 25, 26] suggests that dynamical preservation 
of nonlinear constraint distributions may also be possible. 
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A Proof of theorem 4.3 

We only need to prove the second statement. Since Mx and u-ji do not depend on x, we have 

= (A-.V)i^ + (V.A-)-'' 



Then, upon calculating 

61 „ „ d 51 dU 



5u^ 



O's.dj ax J 
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the nonholonomic part of the first equation becomes (in vector notation) 

d 51 , 6l\ 51 

V £x X cr — / X cr = 

dtdX SX ■' d^6f 



/((X■V)(z/o■(7^)))xo■(7^) + /(^9xf/* j fdwdi^j xa- 
f {{^i.'^n^n-^^)^ + aj^xcr)) xcr + / (^JJ * j f &vdu^ x cr 
fcrxcrxa^ — f(Tx{i/xi/x n(7^)) + / ^y_U * J /dvdi/^ x cr 



where we have used the continuity equation (35). Therefore, the first Euler-Poincare equation 
reads as 

Finally, one computes the micropolar terms (right hand side above) as follows 



A 

T " 



dtJx^^^'Jx)^ -^dnVf^ 

d 61 61 ^ ■ 61 ^^. 61 SI ^. d 61 ^t\^ 

+ ■ ^ + ^ ^ r ^ j ^ + ^ ^^^^ - ^ j 

= f(^j{n)a, + {Vn{i^n)^\jmun)^,)n^ - ^^Tr(z/^j(^)^)^^ - En^+ (d^U * j /dvdi/^ 7^ 

= / (^J(7^)a, - J(7^)^.^. - En^(7^) + (^^Tef/ * j f dvdiv^ 7^^^ 

so that, in vector notation, the first Euler-Poincare equation becomes (36). ■ 
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